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Abstract. We present two versions of third order accurate jet schemes, which achieve high 
£N| order accuracy by tracking derivative information of the solution along characteristic curves. 

t-H For a benchmark linear advection problem, the efficiency of jet schemes is compared with 

WENO and Discontinuous Galerkin methods of the same order. Moreover, the performance 
of various schemes in tracking solution contours is investigated. It is demonstrated that jet 
schemes possess the simplicity and speed of WENO schemes, while showing several of the 
advantages as well as the accuracy of DG methods. 

< 

CM 

<T^ 1. Introduction 

The advection of field quantities under a velocity field is an important sub-problem in 
many computational projects. Examples are the passive transport of concentrations, or the 
movement of interfaces using level set approaches [T3] . We consider the linear advection 
equation 

(f>t + v-V(j) = 0, (1) 

which moves a scalar quantity (ft(x,t) by a given velocity field v(x,t). Equation ([!]) is aug- 
{S| mented with initial conditions (j>(x,0) = <3?(x) and boundary conditions. All data is assumed 

smooth in space and time. Furthermore, it is assumed that the problem at hand can be de- 
scribed by, or embedded into, a rectangular computational domain, equipped with a regular 

grid. 
_j_ ° 

High order accurate numerical approximations of Q on a fixed grid are commonly based on 
schemes that employ polynomials of sufficient degree to interpolate smooth solutions with the 
required accuracy. One popular type of approach are finite difference ENO [TH] or WENO [TT] 
methods. These store approximations of the solution values at the grid points, and achieve a 
high order polynomial approximation by considering local neighborhoods that are several grid 
points wide. Another type of approach are discontinuous Galerkin (DG) [U HI [5] methods. 
These achieve high order accuracy by storing a high degree polynomial approximation in each 
grid cell, and approximating the flux through cell boundaries based on a weak formulation of 
@. Both WENO and DG methods are based on a semi-discretization of Q , and achieve high 
accuracy in time by using strong stability preserving (SSP) Runge-Kutta schemes [TBI [6l 17]. 

Both types of approaches incur problems and difficulties. Due to the wide stencils of 
WENO methods, it is challenging to preserve high order near boundaries. Furthermore, their 
non-locality poses difficulties for an effective parallelization, and their use in conjunction 
with adaptive grid approaches [2j. In contrast, in DG methods, communication is limited 
to neighboring cells. However, DG approaches are characterized by costly quadratures over 



C3 



2000 Mathematics Subject Classification. 65M25; 65M12; 35L04. 

Key words and phrases, jet scheme, advection, high-order, WENO, DG, efficiency, contours. 



1 



2 



P. CHIDYAGWAI, J.-C. NAVE, R. R. ROSALES, AND B. SEIBOLD 



grid cells and their edges, significant time step restrictions, and non-trivial implementation 
aspects. 

A new class of approaches for ([I]) , so called jet schemes, has been proposed recently p2J [15] . 
The goal of these methods is to provide an attractive compromise between WENO and DG 
methods, possessing the optimal locality and good resolution properties of DG, while being 
close to the computational efficiency and ease of implementation of WENO schemes. In this 
paper we conduct a comparative study of the efficiency and accuracy of this new class of 
methods. We consider two versions of third order accurate jet schemes, and compare them 
with WENO and DG methods of the same order. The considered jet schemes are described 
in Sect.[2j Then, in Sect. [3| DG methods are outlined, and in Sect. [I] information about the 
WENO schemes that we employ is given. Sect. [5] shows the numerical results obtained when 
applying the three types of approaches to two versions of a benchmark problem. In Sect. |6j a 
discussion of the observations and theoretical estimates of the computational cost are given. 



2. Jet Schemes 



Jet schemes are based on an advect-and-project approach in function spaces. Given an 
approximation to the solution of ([I]) at time t n , an approximate solution at time t n+ \ = t n +At 
is obtained by the time step 

r +1 = P o A tn+1 , tn r ■ 

Here At n+1} t n is an approximate advection operator, defined by evolving the solution along 
characteristics using a numerical ODE solver, and P is a projection operator, given by a 
piecewise Hermite interpolation based on parts of the jet of the solution at the points of 
a cartesian grid. The fact that the projection requires data only at grid points allows to 
use this approach as a numerical scheme: when evaluating the advection operator, only the 
characteristic curves that go through grid points at t n+ \ need to be considered. In addition, 
the evolution of derivatives of the solution along these characteristic must be found. As 
outlined below, the required spacial derivatives of the advection operator can be found by 
analytical differentiation (Sect. 2.3) or by approximations based on tracking multiple nearby 
characteristics (Sect. 2.5). 

Jet schemes can be constructed in any space dimension, and for any order of accuracy 
|15j . For simplicity, here we only describe third order schemes in two space dimensions. We 
consider a rectangular computational domain Q C M 2 , equipped with a regular cartesian grid 
of grid size h. 



2.1. Projection. In a grid cell [a, a + h] x [b,b + h] C £1, let the vertices be indexed by a 
vector q € {0, l} 2 , such that the vertex of index q is at position xg = (a + hq\,b + hq-z). On 

each of the four vertices, let a vector of data be given by </>| V a G {0, l} 2 . The data represents 
partial derivatives of orders up to 1 in each variable, as follows: a = (0, 0) represents function 
values 4>; a = (1,0) and a = (0,1) represent first derivatives d x (j> and d y 4>, respectively; and 
a = (1, 1) represents d xy 4>. This data is interpolated by the bi-cubic polynomial 

E 4wl(x), (2) 
q,de{o,i}p 

where the Wi{x) are bi-cubic basis functions, given by the tensor product formulas 

Wl(x) = h^ <(^)<(^), 
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and the w% are the univariate basis functions 

w%(x) = 1 - 3x 2 + 2x 3 , Wq(x) = 3x 2 - 2x 3 , w%(x) = x - 2x 2 + x 3 , w\(x) = -x 2 + x 3 . 

The bi-cubic interpolant ^ is an 0(/i 4 ) accurate approximation to any sufficiently smooth 
function (j) that it interpolates on a cell of size h |15j . 

On the computational domain Qcl 2 , let the grid points be labeled x^, where m G Z 2 . 
For any (sufficiently smooth) function : Q — > M, we define a global interpolant as 
follows: at each grid point x^, evaluate the derivatives of 0, d x (f>, d y <fr, and d xy 4>, to produce 
a data vector 0™ Va 6 {0, l} 2 . Then, on each grid cell, use this data to define the bi-cubic 
interpolant ([2]). In the function space 

S 2 ' + = {<;/! G C 1 : if) twice differentiable a.e. with D 2 ip G , 

this procedure can be applied using the following convention: whenever d xy ip must be evalu- 
ated at a point at which D 2 xp is not defined in the classical sense, we define it as 

d X yip(xm) = \ ( ess lim sup d xy ip(x) + ess lim inf d xy ip{x) 



The view on the general order of approximation case [T5] reveals the rationale for the notation 
S 2 ' + . Clearly, the re-application of the interpolation procedure does not change the result: 
T-Lu^ = H<f>- Hence, in the space S 2 > + it can be formulated as a projection operator 

Pcp = n^. (3) 

2.2. Advection. The characteristic form of equation ([!]) is 

d<f> . dx ^ 

— = along — = v(x, t) . (4 
dt B dt K ' y ' 

Let X(x,r,t) denote the solution of the ODE for the characteristic curves at time t, when 
starting with initial conditions x at time r, i.e. it is defined by 

—X(x,T,t)=v(X(x,T,t),t) with X(x, t, t) = X . 

Then due to Q, the solution of satisfies 4>(x, r) = 4>(X(x,T,t),t). In practice, the char- 
acteristic ODE Q must be approximated. Let X(x, r, t) represent an approximate solution 
of the ODE for the characteristic curves at time t, when starting with initial conditions x at 
time r. It typically arises from a numerical ODE solver, e.g. a high order Runge-Kutta step. 
Furthermore, introduce the associated approximate advection operator A T j, which maps the 
solution at time t to an approximate solution at time r. It acts on a function g{x) as follows 

(A Tit g){x)=g(X{x,T,t)) . 

As shown in [12], the use of a locally k th order Runge-Kutta scheme results in a k th order 
accurate approximation to the solution 

A t+A t,t <f>(x, t) - <t>(X(x, t + At, t),t) = 0{\At\ k ) . 

This means that through the characteristic equations Q , any ODE solver induces an approx- 
imate advection operator of the same order of accuracy. However, this idea alone cannot be 
used as a numerical time stepping scheme, since the step from t to t + At in general generates 
a function A^Att 9 that cannot be represented with a finite amount of data. Therefore, at the 
end of every time step, we apply the projection operator ([3]), which generates a function that 
can be stored on a computer. The function space S 2,+ is invariant under diffeomorphisms, 
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thus both P and At+At,t map from S 2,+ into itself. Consequently, one full approximate 
solution step is given by applying P o At+At,t to the approximate solution at time t. 

2.3. Derivative Updates by Analytical Differentiation. The application of the projec- 
tion Q requires the knowledge of ift, d x ip, d y ip, d xy ift at grid points, with ift = A t +At,t <ft n - 
These spacial derivatives of At+At,t 4> n can be found by analytically differentiating the ODE 
solver, as demonstrated below for the Shu-Osher scheme [16] . One step with this scheme is 
0((Ai) 4 ) accurate, which for the scaling At oc h matches the 0(h 4 ) accuracy of the projection 
operator Q. Using the notation t n = t, t n +i = t + At, and t n+ i = t + \At, we obtain the 
updates for the required parts of the jet (ft, S7<ft = (d x <ft, d y cft), and d xy (ft at a grid point x as: 



X\ 


= x - Atv(x, t n+1 ) 


Vxi 


= I — At Vv(x, t n+1 ) 


d xy X\ 


= -Atd xy v(x,t n+ i) 


X 2 


= \x + \x\ - \Atv(x\,t n ) 




= fl+±Wi - iAiVxi ■ W(xi,t n ) 


d xy x 2 


= \d xy xi - \At(d xy xi ■ Vv(xi,t n ) + {{d x x{) T - {d y x\)) : D 2 v(xi,t n )) 


•^foot 


= \x+ §x 2 - lAtv(x 2 ,t n+ i) ^ 


VXfoot 


= |I+ |Vf 2 - | AtX7x 2 ■ Vv(x 2 ,t n+ i) 


d xy X{ QO t 


= ld xy x 2 - lAt(d xy x 2 ■ S7v(x 2 ,t n+ i) + ((d x x 2 ) T ■ (d y x 2 )) : D 2 v(x 2 , t n+ i )) 


(ft(x, t n+ i) 


— %(xf 00 t) t n ) 


(V0)(x, t n+1 ) 


= Vff 00 t • VH(xf oot ,t n ) 


(d xy (ft)(x, t n+1 ) 


= d xy x loot ■ VH(x ioot ,t) + ((d x Xf oot ) T - (dyXfoot)) : D 2 U{x ioot ,t n ) 



In this approach, the characteristic curve is tracked from x at time t + At back to Xf 00 t at 
time t. The data at this position is given by the Hermite interpolation, defined by the data 
(at time t) at the four vertices of the cell that Xf 00 t is contained in. The update rules for the 
derivatives are systematically inherited from the update rule for the function value, and they 
match exactly what one would obtain when evolving the solution using At+At,t everywhere, 
and then applying the spacial derivatives. 

2.4. Order of Accuracy and Stability. One step of the jet scheme described above is 
fourth order accurate, since the advection operator is 0((At) 4 ) accurate, and the projection 
operator is 0(h 4 ) accurate. As usual, when going from this local error to the global error, one 
order is lost, since O(^) time steps are required to reach the final time. Hence, the presented 
scheme is globally third order, given that it is stable. 

Like in many other numerical approaches, stability is not automatically guaranteed. As 
shown in [15], jet schemes can be constructed (by using a different projection than ([3])) that 
are unstable. However, the jet scheme based on the projection Q is stable. A key factor 
in the stability argument is that among all (sufficiently smooth) functions that match given 
data on a cartesian grid, the Hermite interpolant ^ minimizes the stability functional 

F[(ft}= \ (d xxyy (ft(x)) 2 dx . 
Jn 
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Hence J-[Pcj)] < J-[(j>] for all sufficiently smooth (j), which yields bounds on the amounts of 
oscillations that the numerical scheme can create [T3J . 

2.5. A Simpler and More Efficient Derivative Tracking. The analytical differentiation 
procedure, outlined with example ([5]), is the most general approach to derive the update rules 
for derivatives. However, it is not always the simplest to implement. An alternative approach, 
which is solely based on tracking function values, is provided by e-finite differences. For the 
here considered third order jet scheme in 2D, the following procedure can be employed. As 
shall be seen in Sect. [5j it leads to a very simple and efficient numerical scheme. 

At a grid point x = (x,y), instead of tracking one characteristic curve from time t n +i 
back to time t n , we track four characteristic curves, starting at x q = (x + q±e, y + q^e) where 
q £ {—1, l} 2 , and e is a small number, see below. Let the corresponding characteristic foot- 
points be denoted x^ oot . The updates for the required derivatives are obtained as follows: 

(1) The "center of mass" + x^^ + Xfo' ot ^ + ^f 00 t' ^) determines to which cell all 
four foot-points are associated to. 

(2) The Hermite interpolant ([2]) that corresponds to that cell is evaluated at each foot-point 
to yield the values 4>(x q , t n+ i) = H(x^ oot ,t n ). 

(3) The required parts of the jet are defined as 

cf>(x, t n+l ) = \ {4>{x^\ t n+ i) + ^(^^'.Wi) + 4>{^- 1 '- 1 \ t n+1 ) + 4>{x^- l \t n+1 )) 

d x (j>(x,t n+l ) = £ ( < /»(x( 1 ' 1 \t n+ i) - 0(f(- 1 ' 1 ),t n+1 ) + - ^C^" 1 '- 1 ), *„+!)) 

dy<t>{x,t n+l ) = £ ( < /»(x( 1 ' 1 \t n+ l) + 0(f(- 1 ' 1 \t n + 1 ) - 0(#'- 1 ),t n+ l) - ^(^ lrl) ,Wl)) 

d xy <p{x,t n+1 ) = ^{4>{x^ l \t n+1 ) -^(#- 1 - 1 ),t n+ i) -0(f( 1 '~ 1 ),t n+1 ) + ^(f(- 1 - 1 ),t n+ i)) 

All approximations are 0(e 2 ) accurate. In addition, round-off errors of up to 0(5/s 2 ) arise, 
where 5 is the accuracy of the floating point operations. With the optimal choice e = 0(5 1 ^), 
the errors incurred by the e-finite differences are of mag nitude 0(5 1 / 2 ). Thus, with double 
precision arithmetics, the presented approach can be used up to a desired accuracy of 10 -7 . 

3. Discontinuous Galerkin Method 

DG methods |14|, 0] have a wide area of application. They can successfully approximate 
nonlinear problems on unstructured geometries, and thus are much more general than the 
problems considered in this paper. However, since the task at hand is the numerical approx- 
imation of the linear advection equation |T| on a regular grid, it is natural to investigate the 
accuracy and efficiency of DG methods for it. The DG methodology that we use here is in 
line with the ideas presented in [3 [5]. For convenience, we restrict the presentation here to 
incompressible velocity fields V -v = 0. In this case, the advection equation can be written 
in conservative form as 

cj)t + V • (v4>) = . (6) 

Equation ([6| is discretized in space as follows. Let Th be a regular triangulation of the 
computational domain S7. At each instance in time, we seek an approximation 4>h of </>, such 
that 4>h{t) belongs to the finite dimensional space 
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where P$(K) denotes the space of polynomials of degree < k that live on the element K. 
Here, we choose k = 2 to obtain a third order scheme. As in the standard DG formulation 
|3j, we replace (/> by in Q, multiply by a test function ip^ G S*(K), and integrate over 
K £ Th to obtain the semi-discrete formulation 

d 
It 



4> h ip h dx= v-Vip h (f) h dx- (j) h v ■ n e ^ K ip h dT Vip h eSg(K). (7) 

K JK JeedK 



Here n ej x is the outward unit normal on edge e, and 4>h is an expression that comprises the 
values of <ph on each side of the edge. We use an upwind flux, which is defined as 

if v ■ n e ^K > 
otherwise , 

where (fyr = lim e ^ ± <f>h(x + en et K,t), and the velocity field is evaluated in the edge center. 
The integrals over edges and elements in ([7]) are approximated by Gaussian quadrature rules 
that are exact for polynomials of degree 2k over the elements, and exact over polynomials of 
degree 2k + 1 over the edges. This results in an ODE system of the form 

j t 4> h = L h (<i) h ) with A (x,O) = , (8) 

where Lh{4>h) is the spacial discretization of the operator —V • (vcj)), and &h{x) £ Sfc approxi- 
mates 3>(x). The time integration of ^ is done using the third order SSP Shu-Osher scheme 
|16| . Stability is ensured if the specific CFL condition 

At < d2t+T) h 

with c > 1 is guaranteed [5]. For the examples considered in this paper, the choice c = 2 
turns out to yield the lowest cost vs. accuracy ratio. Since for the numerical tests conducted 
here the solutions are very smooth, no slope limiters are implemented. 

The DG code used here is based on triangular elements. In order to have a fair comparison 
with the cartesian grids that WENO and jet schemes use, we design meshes as follows. For 
a desired resolution h, first a cartesian grid of resolution \[2h is constructed. Then each cell 
is divided along its diagonal into two triangles. Like a cartesian grid of resolution h, the 
resulting triangular mesh consists of 1/h 2 elements, and the shortest height in an element is 
h. 



4. WENO Scheme 

Like DG methods, WENO schemes are based on a semi-discretization of ([!]) in space, 
and the use of SSP schemes [16, 6j to advance in time. Unlike DG methods, WENO schemes 
are specific to regular grids. We employ the third order WENO finite difference approach 
described in [9], using the Shu-Osher scheme [16] with At = h in time. For the spacial 
approximation, equation is rewritten as 

4> t = -ud x cf> -vd y (f) . 

At each grid point, both derivatives d x (j) and d y (p are approximated in a univariate fashion 
by using values on four grid points in each coordinate direction: two in the upwind direction, 
and one in the downwind direction. WENO schemes have limiters built into the derivative 
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approximations: smoothness indicators yield a nonlinear weighted average of different polyno- 
mial approximations of the derivatives. We consider two versions of limiting: one as described 
in [9], and another one without limiting (which yields a linear finite difference scheme). 



5. Numerical Results 



We test the accuracy, the relative efficiency, and the quality of tracking contours of jet 
schemes, DG methods, and WENO approaches using the classical vortex in a box flow test 
[UGUI, adapted as follows. On the computational domain (x,y) G [0, l] 2 , and for t G [0,T], 
we consider the linear advection equation ([I]) with the velocity field 

-»/ ja ( t\ ( sin 2 (7nr) sin(27ry)\ fnS 

<*,y,t) = ™ S (nt)^_ sH lJ ) J^ y >), (9) 

which is a model for the passive swirling and successive un-swirling of a concentration field 
by an incompressible fluid motion. This test is a mathematical analog of well-known "un- 
mixing" experiments [SJ [TTJ. For this flow, the final solution equals the initial conditions, 
i.e. 4>(x, y, T) = cp(x, y, 0). We consider smooth initial conditions, and periodic boundary con- 
ditions. Note that for linear advection problems, this test is quite general. Moreover, it has 
a built-in difficulty parameter: the larger T, the more challenging it is to numerically resolve 
the highly elongated contours of the solution. 



5.1. Numerical Efficiency. For the investigation of numerical efficiency, we consider smooth 
periodic initial conditions <f>(x, y, 0) = cos(27rx) cos(47ry), and choose T = 1. For various mesh 
resolutions h, we apply five different numerical schemes to the test problem. The results are 
shown in Figs. fTH3l Specifically, we consider: 



(i) A classical third order WENO scheme (see Sect. [4]), denoted by large squares; 

(ii) A third order WENO scheme without limiting (see Sect. [4]), denoted by small squares; 

(iii) A third order DG method (see Sect. [3]), denoted by triangles; 



(iv) A third order jet scheme with a full derivative update (see Sect. 2.3), denoted by circles; 



(v) A third order jet scheme based on e-finite differences (see Sect. 2.5), denoted by stars. 



For all methods, we measure the error at t = T in the L°° norm, and the CPU time. Figure [T] 
shows the error as a function of the resolution h. One can observe that jet schemes and DG 
have roughly the same accuracy. In contrast, for the same resolution, WENO schemes yield 
significantly larger errors. It is also visible that the WEN03 with limiters does not achieve 
the full third order accuracy on the range of resolutions under consideration (see [9]). It is 
apparent that for the smooth solution at hand, limiters are not beneficial. Figure [2] shows the 
CPU time as a function of the resolution h. One can see a clear ranking in computational 
costs: WENO schemes are very fast, jet schemes are more costly, and DG is even more 
costly. It is visible that e-finite differences are not only simpler to implement, but are also 
significantly faster than the jet scheme based on the full derivative update. Finally, Fig. [3] 
shows the CPU time as a function of the error. This cost vs. accuracy ratio measures the 
true efficiency of the numerical schemes. The results show that the low computational cost 
of WENO schemes renders them preferable over DG. However, WENO does not possess the 
optimal locality that DG has. Jet schemes provide an interesting alternative. While the full 
update version is slightly less efficient than the non-limited WENO, the e-finite differences 
version is even more efficient. 
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Figure 1. Error convergence for third order jet Figure 2. Scaling of the computational cost for 



schemes, DG, and WENO. 



jet schemes, DG, and WENO. 



-O-WEN03 

- ■ -WEN03 no limiters 
-4-DG 

-•- jet scheme full 
f- jet scheme e-FD 




L°° error 



Figure 3. Efficiency of jet schemes, DG, and 
WENO, measured as cost over accuracy. 



5.2. Quality of Contours. In order to assess the quality of the numerical approximations 
that the different methods produce, we consider how accurately contours of the solution are 
deformed under the flow ^ with T = 6. For this test, we choose the initial conditions 

<f>(x, y, 0) = exp (-10(x - 0.5) 2 - 10(x - 0.75) 2 ) , 

and consider three contours at cj) = exp( -Wr 2 k ), where n = 0.044, r 2 = 0.132, and r 3 = 0.220. 
At t = 0, these contours are three concentric circles of radii r^, centered at the point (0.5, 0.75). 
We plot the deformed contours at t = % = 3, the time of their maximum deformation. 
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Figure 4. Deformed contours, computed with Figure 5. Deformed contours, computed with 
a jet scheme using h — 1/90. a DG scheme using h = 1/90. 




Figure 6. Deformed contours, computed with 
a WENO scheme using h = 1/90. 



We compare a jet scheme, DG, and WENO for the same resolution h. As found in Sect. 5.1 



these methods involve very different computational costs. Yet, this equal-/i comparison is of 
practical relevance, since frequently the advection equation ([I]) is only a sub-problem in a 
larger project, and the computational grid is determined by the underlying application. Fig. 
shows the contours obtained with a jet scheme based on e-finite differences (see Sect. |2.5 
Fig. [I] shows the contours obtained with DG (see Sect. [3]); and Fig. [6] shows the contours 
obtained with WENO without limiting (see Sect. [4]). In all cases, the inner (red) contour is 
4> = exp(— Wrf), the middle (purple) contour is <fi = exp(— 10rf), and the outer (blue) contour 
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is cj> = exp(— 10 r|). Moreover, each contour underlays the contour of the true solution as a 
thick gray curve. In addition, the computational grid is shown in the background. Note 
that in all cases, contours are plotted using a sub-grid: for jet schemes, the solution inside 
a grid cell is given by the Hermite bi-cubic interpolant ([2]); for DG, the solution inside an 
element is given by a bivariate quadratic polynomial; and for WENO, we use a simple bi- linear 
interpolant on each grid cell. 

As it was already visible in the error convergence graph shown in Fig. [TJ also in terms of the 
tracking of contours, the jet scheme as well as DG are significantly more accurate than WENO. 
The comparison between jet schemes and DG is more interesting. Clearly, the DG contours 
shown in Fig. [5] are closer to the true solution's contours than the jet scheme's contours are. 
However, at the tips of contours, DG produces small ripples and break-ups. In contrast, 
the jet scheme contours are extremely uniform, all the way to their end. Moreover, for the 
tracking of contours, DG has another potential disadvantage: due to the discontinuities of the 
numerical solution across element boundaries, the resulting contours may have small jumps 
(in fact, the contours in Fig. [5] possess jumps, but they are too small to be visible to the eye). 
In certain applications, such as level set methods [13] , this could lead to problems. These 
observations come in addition to the fact that for the same resolution, DG is significantly 
more costly than jet schemes (see Fig. [2]). 



6. Discussion 

The two types of traditional numerical methods, WENO and DG, achieve high order ac- 
curacy in very different ways: WENO considers function values in a wide neighborhood. In 
contrast, DG uses high order polynomials in each element. This gives DG a certain level of 
sub-grid resolution, and restricts communication to neighboring elements only. While based 
on a very different methodology, the recently proposed \12\ [T5] jet schemes share these ad- 
vantages with DG. 

The results presented here show that the advantages of DG come at the expense of efficiency. 
This observation can be explained by the following theoretical estimates. With a third order 
Runge-Kutta scheme, all presented schemes require three right hand side evaluations per time 
step. On a mesh of N elements, even a fully optimized and problem-specific DG code requires 
at least 11.5iV evaluations of the velocity field per right hand side evaluation (7 quadrature 
points on each of N elements, and 3 quadrature points on each of |iV edges). In comparison, 
WENO required N velocity field evaluations, and a jet scheme (using e-finite differences) 
requires AN velocity field evaluations (4 characteristics per grid point). In addition, due to 
its more restrictive CFL condition, DG requires 5-10 times more time steps than WENO or 
jet schemes. The results shown in Fig. [2] are in line with these theoretical estimates, albeit 
at more extreme ratios between the costs of the methods. These increased ratios are most 
likely due to additional data structure management requirements for DG, and also for jet 
schemes. However, the key results in the efficiency plot shown in Fig. [3] remain valid even if 
the measured CPU times were replaced by the theoretical cost estimates given above. 

In terms of accuracy, jet schemes and DG yield similar accuracies (for the same number of 
elements). In contrast, the low cost of WENO is outweighed by its lower accuracy. Therefore, 
in terms of true efficiency, jet schemes and WENO are in the same range, while DG methods 
are significantly more costly. This observation is further confirmed by a comparison of the 
quality of contour approximation. 
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These investigations lead to the conclusion that jet schemes indeed fill a need as compu- 
tational approaches for advection problems on regular grids, namely: they are simpler to 
implement, and less costly than DG. However, in contrast to WENO, they are optimally 
local and thus are preferable for the purpose of parallelization, adaptive mesh refinement, 
and the implementation of complicated boundary conditions. In addition, jet schemes impose 
no restrictions on the ODE solvers that are needed. In contrast, DG and WENO methods 
require SSP ODE solvers to ensure TVD stability [Tj- The key weakness of jet schemes is their 
limited range of applicability. While established for linear advection |121I15|. their application 
to more general advection problems is a the topic of current research. 
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